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Abstract 

Context. Cosmic microwave background (CMB) data analysis 

Aims. Develop and validate tools to estimate residual noise covariance in Planck frequency maps. Quantify signal error effects and 
compare different techniques to produce low-resolution maps. 

Methods. We derive analytical estimates of covariance of the residual noise contained in low-resolution maps produced using a 
number of map-making approaches. We test these analytical predictions using Monte Carlo simulations and their impact on angular 
power spectrum estimation. We use simulations to quantify the level of signal errors incurred in different resolution downgrading 
schemes considered in this work. 

Results. We find an excellent agreement between the optimal residual noise covariance matrices and Monte Carlo noise maps. 
For destriping map-makers, the extent of agreement is dictated by the knee frequency of the correlated noise component and the 
chosen baseline offset length. The significance of signal striping is shown to be insignificant when properly dealt with. In map 
resolution downgrading, we find that a carefully selected window function is required to reduce aliasing to the sub-percent level at 
multipoles, £ > 2A^ sid e, where A^de is the HEALPix resolution parameter. We show that sufficient characterization of the residual noise 
is unavoidable if one is to draw reliable contraints on large scale anisotropy. 

Conclusions. We have described how to compute the low-resolution maps, with a controlled sky signal level, and a reliable estimate 
of covariance of the residual noise. We have also presented a method to smooth the residual noise covariance matrices to describe the 
noise correlations in smoothed, bandwidth limited maps. 

Key words, cosmic microwave background — Cosmology: observations — Methods:data analysis — Methods: numerical 



1. Introduction 



Over the last two decades observations of the cosmic microwave background (CMB) have led the way towards the high precision 
cosmology of today — a process best emphasized recently by the high quality data set delivered by the WMAP satellite. The next 
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major and nearly imminent step in a continuing exploitation of the CMB observable will be data analysis of data sets anticipated 
from another satellite mission, Planck. Planck will observe the entire sky in multiple frequency channels, promising to improve 
over the recent WMAP constraints on many fronts. In particular Planck, as a satellite, will provide us with a unique access to the 
largest angular scales, in which the total intensity has proven controversial and difficult for theoretical interpretation and is still 
poorly measured and exploited in the polarization. Planck will be the only CMB satellite deployed in the next decade. It is therefore 
particularly important that the constraints at large angular scale derived from the anticipated Planck data are not only robust but 
also efficiently exploit the informat ion contained in them This will be certainly necessary if Planck is to set strong constraints on 
the CMB B-mode power spectrum (lEfstathiou et al.l 2009) — one of the most attractive potential science targets of the mission. 

The analysis of constraints on the largest angular scales requires robust statistical estimators accounting for a proper description 
of the statistical properties of the sought-after sky signal, instrumental noise and other residuals due to the instrument, astrophysical 
signals and/or data processing. This paper focuses on two of those ingredients — instrumental noise and so-called pixel noise, the 
latter due to residual sky power on sub-pixel scales. In the standard data analysis pipeline the measured time ordered data are first 
projected onto the sky, an operation called map-making, producing map-like estimates of the sky signal, which are subsequently 
analyzed, e.g., in order to derive constraints upon the power spectrum. The map-making process is usually understood to be a linear 
operation on the input, measured data and therefore a statistical uncertainty of the produced sky maps can be straightforwardly 
obtained given known characteristics of the time-domain data. Those also usually involve assumptions about piece- wise stationarity 
of the instrumental noise, assumed to conform with Gaussian statistics. In the realm of the Planck analysis such a straightforward 
route is not however plausible. This is because of the high resolution of the Planck instruments, the full sky coverage and length of 
the mission combined with the high sampling rate of the sky signals. That leads to the data set which is large in terms of the numbers 
of both the sky pixel s contained in the maps and the directly regi st ered m easurements. The map-making procedures developed in the 
context of Planck (iPoutanen et al.ll200q: lAshdown et al . 2007a. b. [20091) have been demonstrated to be capable of dealing with the 
expected volumes of the data, producing high-quality maps; nevertheless the calculation or even just storage of their respective noise 
covariance matrices at their full resolution is beyond the limits of even the largest currently available supercomputing facilities. This 
is because, unlike maps — sizes of which scale linearly with a number of pixels, Af p i x , — the noise covariance matrices scale as a 
square of it and their inversion involves 0(N^ ix ) floating point operations. We emphasize that due to a combination of its scanning 
strategy and noise-like contributions correlated over long periods of time, we expect that non-negligible large scale noise correlations 
will be present in the maps derived from the Planck data and will be particularly important in the analysis of the polarized signals 
given their lower amplitudes. 

In this paper we develop tools necessary for the statistically sound analysis of constraints on large angular scales. These include 
robust approaches to producing low-resolution maps and techniques for estimating pixel-pixel correlations due to their residual 
noise. The low-resolution maps are expected to compress nearly all the information relevant to the large angular scale in fewer 
pixels and are therefore more readily manageable. Given our Gaussian noise assumption the full statistical description of the map 
uncertainty is given by its pixel-pixel noise covariance matrix (NCM). This is defined as 



and s is the 3A^ p i x input sky map of Stokes I, Q, and U parameters and s is our estimate of s. In the absence of signal errors, the 
difference, (s - s), contains only instrument noise. We note that N is a symmetric and usually dense matrix, which in general will 
depend on the map-making method that produced the estimate. In the following we will consider a number of numerical approaches 
to calculate such a matrix for each of the studied low-resolution maps and then test their consistency with the actual noise in the 
derived maps. 

The full noise covariance matrices have been commonly computed and used in the analysis of the sm all-scale balloon-borne, 
(e.g. lHananv et alJl2000b Ide Bernardis et alJl2000h and ground-based experiments (e.g ., iKuo et al.ll2004l) . The COBE-DMR team 
also used them to derive low-^ constraints, ( e.g.. lG6rski et al.lll996UWright et al.lll996l) . In all those cases, however, no resolution 
downgrading has been required, unlike with Planck, as the calculations for those experiments could be done at a full resolution. 
To this date, only the WMAP team has encountered a similar problem. The instrument noise model employed by t hem is in fact 
similar to t he one used for Planck. It is parametrized, however, in the time domain rather than in the frequency domain (iJarosik et al.l 
l2(Ml2007h . Calculation of the WMAP NCM is formulated in exactly the same manner as for our optimal map-making method 
and the WMAP likelihood codeQ ships with an NCM very similar to what we present here, although without the II, IQ and IU 
covariance blocks. The simplification is motivated by the high S/N ratio of the low temperature multipoles and weak coupling 
between temperature and polarization pixel noise. 

Our analysis is made unique by the differences in the experiment design: WMAP pseudo-correlation receivers are differencing 
assemblies (DA) with two mirrors, whereas Planck will use a single mirror design (HFI, the high frequency in strument) or has 
a reference load in place of the second mirror (LFI, the low frequency instrument) (iPlanck Collaboration] 120051) . Between these, 
the pixel-pixel correlations are different. In principle the balanced load systems of COBE and WMAP should bring less correlated 
noise than the unbalanced Planck LFI. On the other hand, differencing experiments generate pixel-pixel noise correlations even 
from white noise, whereas in Planck they originate from the correlated noise alone. In add ition we also study here the so-called 
destriping algorithms, which have be en proposed as a Planck- specific map-making approach (Delabrouille 1998; Maino et al . 2002; 
iKeihanen et aDl2004t iKeihanen et al.ll2005h . 

Residual noise covariance for PL ANCK-like scanning and instrument n oise has been studied b efore, either via some simplified 
toy-models (iStompor & Whitel 12004) or in more realistic circumstances (Efstathiou 2005, 2006). Those studies approached the 
problem in a semi-analytic way and thus needed to employ some simplifications, which we avoid in our work. They were also based 




where ((s - s)) = 0, 



(1) 



1 http://lambda.gsfc.nasa.gov/ 
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solely on the destriping approximation, assuming that noise can be accurately modeled by relatively long (one hour) baseline offsets 
and white noise, and did not consider any other approaches. In this work we extend those analyses into cases where modeling the 
noise correlations requires shorter baselines and compare those with optimal solutions. 

2. Algebraic background 

2.1. Maps and their covariances 

To formulate map-making as a maximum likelihood problem we start with a model of the timeline: 

d = As + Bx + n, (2) 

where the underlying microwave sky signal, s, is to be estimated. Here A is the pointing matrix, which encodes how the sky is 
scanned and n is a Gaussian, zero mean noise vector, x denotes some extra instrumental effects, usually taken hereafter to be 
constant baseline offsets, which we will use to model the correlated part of the instrumental noise and B — a 'pointing' matrix for x 
describing how it is added to the time domain data. Convolution with an instrumental beam, assumed here to be axially symmetric, 
is already included in s. 

The signal part of the uncalibrated data vector, d, is the detector response to the sky emission observed in the direction of pixel 
p. For a total power detectors, e.g., the ones on Planck, it is a linear combination of polarized and unpolarized contributions: 

d t = 7C[(l + e)s p i + (1 - e) (s pQ cos(2^) + s p u sin(2^))} + n u (3) 

where it is implied that sample t is measured in pixel p andxt is a detector polarization angle with respect to the polarization basis 
and we have dropped the baseline term for simplicity. Eq. [3] includes an overall calibration factor *K, and a cross polar leakage 
factor 6, however, in what follows, we only consider the case of perfect calibration, setting *7C = 1 with no loss of generality, and no 
leakage, 6 = 0. 

To simplify future considerations we introduce a generalized pointing matrix, A', and a generalized map, s' . They are defined 

as 



A' = [A,B] , and s f = 



(4) 



Using those we can rewrite our data model in a more common form, 

d = A V + n. (5) 

The detector noise has a time-domain covariance matrix N = (nn T ) and the probability for the observed timeline, d, becomes the 
Gaussian probability of a noise realization n = d - A's'\ 



P(d) = P{n) = [(27r) dimw det n] m exp f-^AT 1 w j P (x) , 



(6) 



where the last factor is a prior constraint on the noise offsets, x, which hereafter we will take to be a Gaussian with a zero mean and 
some correlation matrix, *P, i.e., 

P (x) oc exp l-ijt 1 ^" 1 *! (7) 

By maximizing this likelihood with respect to the sky signal and baselines contained in s', we find an expression for a maximum 
likelihood estimate which reads, 

s' = (ft" 1 + A /T A/ r_1 A , )~ 1 A /T 7V _1 J, (8) 
where ft -1 is defined as, 





P~ l 



(9) 



The first part of the vector s' is an estimate of the actual sky signal, s, while all the rest is an estimate of the baseline offsets, Jc. The 
map A' T N~ l d is called the noise- weighted map. In case of flat prior for x, expressions identical to Eq. (8]) can be also derived from 
minimum variance or generalized least square considerations and we will refer to s as either a minimum variance or optimal map 
in the future. We note that we have ignored here any pixelization effects that cause differences between d and A's' even for a noise 
free experiment. This is usually true in the limit of the pixel size significantly smaller than the beam resolution of the instrument. If 
this condition is not fulfilled, the pixelization effects may be important and special methods may be needed to minimize them. We 
discuss specific proposals in Sect. 12.31 In the absence of such effects, the difference between a map estimate, and the input map, 
s\ is called residual pixel domain noise. 
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Let us now consider first the pref actor matrix in Eq. ([8]), 
W = + A ,T 7V _1 A r ) _1 , (10) 

a weight matrix combining both the baseline prior and the noise variance weights. It acts on the generalized noise-weighted map, 
producing estimates of the pixels and baselines. 

Given that our generalized map is made of two parts: the actual sky signal and the baseline offsets, the matrix M f has four 
blocks: two diagonal blocks, denoted M and M x , and two off-diagonal blocks, each of which is a transposed version of the other 
and one of which is referred here to as M . Using inversion by partition we can write an explicit expression for each of these blocks. 
For example, for the sky-sky diagonal blocks we get, 



M 



A T N~ l A - (a t N~ 1 b) (<P~ l + B T N l B) 1 (A T N~ l Bf 



(11) 

= (a t N~ 1 A)"' + (A T 7V _1 A)"' (a t N~ 1 b) M x (A T N~ l By (a t N~ 1 a)" 1 , (12) 
while for the offset-offset part, 

P~ l + B T N l B - (b t N~ 1 a) (A T 7V _1 A)" 1 (b t N~ 1 A) 



(13) 



With help of these equations we can now write explicit separate expressions for the estimated sky signal and offsets. The former is 
given by, 

s = (MA T + M B T ) N~ l d (14) 
while the latter, 

x = (M T A T + M X B T ) N~ l d. (15) 

We can also combine these two equations to derive an alternative expression for the sky signal estimate, which makes a direct use 
of the offsets assumed to be estimated earlier, 

s = (A 7 //" 1 A)" 1 (a t N~ 1 d - A T N~ l Bx) (16) 

If the assumed data model, Eq. d2]), and the time domain noise and baseline covariances are all correct, then the covariance of 
the residual pixel domain noise is 

N' = <(S' - s')(? - s'f) = M'. (17) 

In particular, Eq. (Q]), the pixel-pixel residual noise covariance matrix, N, is equal to M and given by Eq. (TT2b . 

We note that a sufficiently high quality estimate of the inverse time domain correlations, TV -1 , is required in order to calculate 
both the minimum- variance map and its noise covariance. If it is misestimated the map estimate will still be unbiased, though not 
any more minimum variance or maximum likelihood, and its covariance will not be given any more by Eq. (fTTl) . 

For example for computational reasons we will find later that using some other matrix, denoted here as At -1 , rather than the 
actual inverse noise covariance, TV" 1 , in the calcul ation of the map e stimates in Eq. (8]) can be helpful. The corresponding noise 
correlation matrix for such a map is then given by (IStompor et al .1120021 ungeneralized case) 

N' = <(S' - s')(s' - s') T > = (ftT 1 + A' T M~ l A'y l (ftr l ft N ftr l + A ,T M~ l NM~ l A') (ft~ l + A' T M~ l A')" 1 , (18) 

where M and ft define our map-making operator, whereas N and ft^ are the true noise properties. This expression is significantly 
more complex and computationally involved than Eq. (fTTl) . Fortunately, as we discuss in the following, in many cases of interest, 
the latter expression turns out to be a sufficiently good approximation of the former with N~ l replaced by M~ l at least for some of 
the potential applications. 

The Planck Worki n g Grou p 3 (CTP) has performed extensive studies of different map-making approaches ( Po utanen et al.1 2006; 
lAshdown et al. 2007a b,|2009|). They have been shown to produce different residual noise structures in the computed maps studied 
in detail in those papers. A map-making method should only be considered complete once the residual noise covariance associated 
with it can be understood and sufficiently well characterized. 

The map-making methods considered for Planck fall into two general classes both of which are described by the equations 
derived above. The first class, called optimal methods, assumes the sufficient knowledge of the time domain noise correlations, does 
not introduce any extra degrees of freedom, x, (or alternately assumes a prior for them with V vanishing). In this case one can derive 
from Eqs. (fT6l) & (fl2l). 



s = (a t 7V _1 a) A T N~ l d, (19) 
N = (A T N- l A)~ l (20) 
where the latter is a covariance of the residual noise of the former. 
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The second class of methods introduces a number of baseline offsets with a Gaussian correlated (so-called generalized destripers) 
or uncorrelated (standard destripers) prior on them and describe the noise as an uncorrelated Gaussian process. The destriper maps 
are evaluated via Eq. ([14l) or Eq. ([T6b , with N assumed to be diagonal. Clearly on the time-domain level the destriper model is just 
an approximation, therefore at least a priori we should use the full expression in Eq. (IT8l) to estimate its covariance. The CTP papers 
have shown that for Planck the two methods produce maps which are very close to one another. Moreover, they have shown that 
using a generalized destriper, the derived maps eventually become nearly identical to those obtained with the optimal methods, if an 
appropriate length of the baseline and a number of the baseline offsets is adopted together with a consistently evaluated prior. This 
motivates using the simplified Eq. (Tf2l) as an approximation for the noise covariance of the destriped maps, i.e., N = M. We will 
investigate the quality and applicability of this approximation later in this paper. 

In this paper we extend the analysis presented in those earlier CTP papers. We first study the covariances derived for the different 
map-making algorithms using Eq. (TT2lh compare their properties and test how well they describe the residual noise in the actual 
maps. As all those calculations can not be performed at the full instrumental resolution, we also discuss methods of producing the 
low-resolution version of the maps. 



2.2. Time domain noise 

We assume that the time domain noise is a Gaussian process and for the simulations we take the noise power spectral density to 
have the form 

2 fa + fa 

P(f) = — — • Jmm J , (21) 

J sample J knee t J 

where the shape is defined by the slope, minimum and knee frequencies (a, / mm and /knee respectively) and the scaling by the 
white-noise sample variance and sampling frequency (<x and /sample)- Two examples of the theoretical and simulated noise spectra 
can be seen in Fig.Q] 

In the calculations of the maps using the optimal algorithms or generalized destripers we will assume that noise power spectrum 
is known precisely. As the noise simulated in the cases analyzed here is piece- wise stationary, with no correlations allowed between 
the data in the different pieces (see Sect. O the respective noise correlation matrix, N, is block Toeplitz with each of the blocks, 
describing the noise correlations of one of the stationary pieces, defined by the noise power spectrum. Given that we will approximate 
the inverse of N as also a block Toeplitz matrix with each blocks given by an inverse noise power spectrum. Though this is just an 
approximation it has been de monstrated in the past that it performs exquisitely well in particularly in the cases with long continuous 
pieces of the stationary noise (Stomp or et al.l l2002). as it is a case in all simulations considered here. 



2.3. Low-resolution maps 

Planck will produce maps with resolution of ~ 5 arc minutes at frequencies of 217 GHz and above, and < 13 arc minutes from 
70-143 GHz. The sky maps pixelized at the full available resolution will therefore include as many as O(l0 7 ^ pixels per Stokes 
parameter. Though it has been demonstrated in previous CTP papers that a calculation of such maps is feasible, the computations 
of the covariances of such maps is clearly well beyond the reach of the current and near- future supercomputers. At the same time 
production of low-resolution maps from data of a high-resolution experiment is not a straightforward task, which in the CMB 
context is made even more difficult due to a disparity in amplitudes of the total intensity anisotropies on the one hand and Q and 
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U Stokes parameters (or E and B polarization modes) on the other. As we emphasized earlier the map-making methods described 
in the previous sections work very well but only in the limit of sufficiently small pixels. Those need to be much smaller than the 
typical variability scale of the considered sky signal, which is usually set by the experimental resolution. Such an assumption is 
clearly not fulfilled in the case of the low-resolution map-making. We therefore expect that the pixel effects are non-negligible in 
this latter case. Moreover as we solve simultaneously for all three Stokes parameters even relatively mild pixel effects present in the 
total intensity maps may have significant consequences for the Q and U Stokes parameter maps. 

In this Section we define three alternative methods of producing low-resolution maps from high-r esolution observatio ns. The 
first two, direct and (inverse) noise weighting methods, have already been used in the WMAP analysis (iJarosik et al.ll2007l) . As the 
third option we consider at the end smoothed (low-pass-filtered) maps and their noise covariance. 

2.3.1. Direct method 

The most straightforward method to produce a low-resolution map is to project the detector observations directly to the pixels of the 
final target resolution. Hereafter we will refer to it as the direct method. 

The direct method is clearly the best choice as far as the described noise covariance is concerned. However, it does not pay any 
particular attention to minimizing the pixel effects. In particular, it may lead to a position-dependent signal smoothing due to a non- 
uniform sampling of the low-resolution pixels — an effect which may further cause aliasing problems at the, for example, power 
spectrum estimation stage. Moreover, for the destripers the direct method means that the baseline offsets are solved at the low target 
resolution. If the subpixel structure of the pixels can be neglected, this will lead to a better determination of the baseline offsets, 
and less residual noise, as the number of crossing points between baselines increases ( Ashdo wn et al.ll2007ab . If, however, sub-pixel 
power is present, it may affect adversely the offset estimation, with magnitude of the effect increasing with the pixel resolution. 
None of the discussed map-making methods is designed to correct for subpixel structure. Therefore the direct method can be taken 
to regard the sky as already smoothed to eliminate the subpixel structure within the large, low-resolution pixels. 

In Sect. [63] we quantify the signal error for the different low-resolution maps and map-makers. 

The noise covariance matrices for such low-resolution maps can be computed directly using formalism presented in Sect. I2.ll 
for example, Eqs. dT71) and dTSb . 

Hereafter, we will use the direct method as a reference with respect to which we compare the other approaches. 

2.3.2. Inverse noise weighting (INW) 

In the case of nested pixelization schemes, such as HEALPix dGorski et al.ll2005l) used in this paper, to downgrade a temperature 
only map, one may compute a weighted average of the subpixel temperatures. A natural choice are the optimal weights, where 
the temperature of a small pixel is weighted with the inverse of its noise variance (or with its hit count provided that the detectors 
have equal noise equivalent temperatures). This weighting leads to the lowest noise of the large pixel in the absence of pixel-pixel 
correlations. We will refer to these maps as inverse noise weighted (INW) maps. 

A similar weighting scheme exists for the polarized data as well. The procedure goes as follows: first the estimated high- 
resolution maps are noise- weighted, then their resolution is downgraded, and the resulting low-resolution, noise- weighted maps are 
subsequently multiplied by the low-resolution noise covariance, which needs to be estimated in parallel. Algebraically, the entire 
procedure can be summarized succinctly on the map level as, 



where W and W" are weight matrices for the high and low resolution maps, respectively. They depend on A and A" — the pointing 
matrices at high and low resolution. X simply sums the pixels in resolution N p [ x to resolution A^ ix , 

Y _ j 1, P subpixel of q m . 
Aqp ~ \ 0, otherwise (ZJj 

In the following we will assume either block-diagonal or diagonal weighting. In the former case the weights are given by, 



while in the latter case they are made of the diagonal elements of the above matrices. Matrix N u is the time domain covariance 
matrix of the uncorrelated part of the noise, n. In the block-diagonal case the noise weighting mixes different Stokes parameters, 
while in the diagonal one each Stokes parameter map resolution is downgraded independently. Throughout this work we will use 
only the block diagonal weighting which, in the cases studied here, turns out to be very close to the diagonal one. 

The covariances for the maps obtained via such a procedure can be derived from Eq. (l22l) and the expressions described in 
Sect.EQ 

For the destriper technique there is one more extra factor which makes this manner of resolution downgrading differ from the 
direct method of the previous Section. As the maps outputted directly by a destriper code are of a high resolution, the baseline offsets 
are also determined at that resolution. If the block-diagonal weighting is then used to downgrade the map, the result is equivalent to 
the direct calculation of the low-resolution map with the baselines determined from the high-resolution analysis, Eq. (TT6b . 



s f = W'XWs, 



(22) 




(24) 



(25) 
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Noise weighting reduces signal errors by first solving the map at a resolution where subpixel structure is weak. In comparison 
to the direct method, the more accurate signal is gained at the cost of higher noise. Like the direct method, INW also disregards any 
subpixel structure but at high resolution. In this case the instrument beam naturally smoothes out small scale structure causing the 
approximation to hold. 



2.3.3. Harmonic smoothing 

Both methods described above may result in the signal smoothing (or its band-width) being position-dependent as it is achieved 
via averaging of the observed high-resolution pixel amplitudes contained within each low-resolution pixel. This may result in the 
aliasing of sky power. 

Applying a smoothing operator to each of the high-resolution maps prior to resolution downgrading could alleviate such a 
problem. The smoothing operation needs to take care properly of the high frequency power contained in the maps avoiding thus 
its being aliased to the power at the scales of interest. As the smoothing operation is usually performed in the harmonic domain it 
requires that the high-resolution map is first expanded in spherical harmonics. If the map has unobserved pixels, they will induce 
undesired mode coupling. For sufficiently complete sky coverage we can "patch" the high-resolution map by adding averages of 
the neighbouring pixels into blank pixels. If the coverage is more incomplete, the missing pixels can be replac ed by a constraine d 
realization of signal and noise, such methods are used for example in the so-called sampling techniques (e.g. J Jewell et al] 12004). 
which have been successfully applied to simulated Planck data. Simple patching will clearly affect only very small scale statistical 
properties. If a constrained realization is applied, the spherical expansion will depend on the input model. In this work we only deal 
with a complete sky coverage leaving an investigation of those effects to the future work. 

To suppress small angular scale power the ex pansion is convolved with an axially symmetric window function (e.g. a symmetric 
Gaussian window function (Chal linor et al.ll200Ql) ). 

aL = W^L W e = e-y^ 2 (26) 
at = 2W ( al, 2 W ( = e~^W, (27) 

chosen to leave only negligible power at angular scales that are not supported by the low target resolution. Finally the regularized 
expansion is synthesized into a low-resolution map by sampling the expansion values at pixel centers. We conduct most of our 
studies using a beam having a full width at half maximum (FWHM) of twice the average pixel side. For the A^ide = 32 resolution 
this is approximately 220' (3°. 7). Whenever transforms between harmonic and pixel space are conducted, it is important to consider 
the range of multipoles included in the transformation. We advocate using such an aggressive smoothing that the harmonic expansion 
has negligible power beyond £ = 3N S id Q and results are stable for any £ max beyond this. For completeness we have set ^ max = 4A^ S id e 
but stress that any residual power beyond £ = 3N S id e will lead to aliasing. 

The sm oothing window does not need to be a Gaussian but it is preferable to avoid sharp cut-offs that may induce "ringing" 
phenomena. iBenabed et al.l 12009) suggest a window function that preserves the signal basically unchanged until a chosen threshold 
and then smoothly kills all power quickly above that angular resolution. Their window is 



1, l<h 

_ [ 1 _1_ ™e ((P - P.\>rrl(P~ - P *\W P. ^ 



Wi = {±[l + cos ((£ - " A))] , h < £ < h (28) 

0, £ > h 

with the typical choice £\ - 5N S i^ Q /2 and £2 = 3N s ^q. 

This method can be considered optimal from the (large-scale) signal viewpoint; however it may be suboptimal as far as the noise 
is concerned, in particular in cases with a strongly inhomogeneous noise distribution on the observed sky. The noise covariance 
matrices described in Sec. l2.1l need to be amended to accurately characterize the residual noise of the smoothed maps and thus we 
need to smooth the matrices as well. 

Smoothing of a map is a linear operation. For any linear operator, L, acting on a map, m, we can compute its covariance as 

{(Lm)(Lm) T ) = L(mm T )L T = LNL T = ^ X { Le t • (L^) T , (29) 

i 

where X{ and e\ stand for eigenvalues and eigenvectors of the noise covariance, N, i.e., N = Yji^fifi\ \ an d m is understood here 
to contain the noise only. We note that in general one should replicate the same processing steps as are to be applied to maps 
and therefore the smoothing operation should be applied to the noise covariance of the high-resolution map and its resolution 
downgraded later. All these steps are described by the operator L introduced above. In this case L is a rectangular matrix with many 
fewer rows (given by the number of low-resolution pixels) than columns (the number of high-resolution pixels). In such a case 
rather than performing the eigenvalue decomposition as suggested by the right-most term of Eq. ([29b , which would require as many 
operation as the cube of the number of high-resolution pixels, it may be more efficient to perform the matrix-matrix multiplications 
in Eq. d29b explicitly. In fact in the latter approach one could rephrase the problem as a series of PCG solutions of a map-making 
type, each of which would result in the computation of the high-resolution covariance, N, times one of the columns of the smoothing 
operator, L T . This could bring the cost of the covariance smoothing down to that comparable with actual map-making operation 
repeated for each of the low-resolution pixels. Though this may be in a realm of capabilities of the present day supercomputers it is 
certainly a huge effort not warranted at the present stage of this investigation. 

Alternately, one may choose to commute the order of the smoothing and downgrading operations as highlighted above. Though 
these two operations are clearly not exchangeable on the map level, due to the potential presence of the sub-pixel power, such an 
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approach can be more justifiable for the noise co variances. In this case we can explicitly compute the low-resolution unsmoothed 
covariance matrices directly and subsequently smooth them with the signal smoothing kernel. One side advantage of this approach 
is that the low-resolution maps are more likely to be genuine full-sky maps than their high-resolution counterparts. Applying the 
smoothing is therefore less likely to require any additional pre-processing. 

In the following we will apply the smoothing technique to both high and low resolution maps already downgraded using some 
of the other approaches. We will demonstrate that such a combined approach results in controllable properties of the residual noise 
on the one hand and well defined sky signal bandwidth on the other. Unlike both the direct and INW methods, the low-resolution 
maps are actually solved from a signal that lacks subpixel structure. 

3. Numerical calculations of residual noise covariance 

This section presents numerical methods to compute the residual noise covariance matrix and describes briefly their implementa- 
tions corresponding to three different map-making methods, the optimal method (MADping and ROMA implementations) and the 
generalized (Madam) and classical destriping (Springtide) methods. 

3.1. Optimal map covariance 

The noise covariance for the maps computed by optimal algorithms using true time domain correlations is given by Eq. (l20l) . The 
calculation of such a matrix proceeds in two steps and two different implementations have been developed in the course of the 
work described here. During the first step the inverse covariance matrix, A T N~ l A needs to be assembled and subsequently inverted. 
Given that the matrix can be singular the latter step needs to be taken with care and a pseudo-inverse may need to be computed. The 
computation of the latter involves a eigenvalue decomposition of the inverse noise matrix. Because the noise matrix is symmetric 
and in principle non-negative definite, its eigenvalues are real and non-negative (^ > 0), and its eigenvectors form a complete 
orthogonal basis. This allows us to expand the matrix as 



Here X\ are the eigenvalues and e\ are the corresponding eigenvectors of the matrix N . We can now invert N 1 by using its 
eigenvalue decomposition, 



and controlling the ill-conditioned eigenmodes. Any ill-conditioned eigenmode will have an eigenvalue several orders of magnitude 
smaller than the largest eigenvalue. By including in the sum only the well-conditioned eigenmodes we effectively project out the 
correlation patterns that our methods cannot discern. This way of calculating the noise covariance is implemented in the ROMA and 
MADping codes. 

MADping is one of the codes of the MADCAF0 s uite of CMB analysis tools. The code is parallelized using MPI and all the 
operations are distributed across multiple processors dBorrilll [T999b . It uses the M3 library mentioned in (Cantal upo et al.ll2QQ9l) 
for data reading and time-domain noise correlation generation. Load balancing is performed based on both the number of pixels 
per processor and the number of time samples falling into those pixels. Each processor scans through its sections of time-ordered 
data (correctly handling overlap with other processors' data) and accumulates its local piece of the inverse noise covariance. These 
pieces are then gathered and written to disk. The scaling of this technique is 



where ^correlation is the filter length set by the noise autocorrelation length. 

For a Planck- sized, full- sky dataset and using a reasonable pixel resolution (half a degree), the construction of the inverse 
noise covariance dominates over the computational cost of inverting this matrix. Nevertheless, inversion methods as well as the 
eigendecomposition scale as 0(N^ ix ). 

In order to correctly treat the signal component of the data map in Eq. ([8]), we must apply our low-resolution noise covariance to 
a noise- weighted map which has been downgraded from higher resolution. This downgrading process is equivalent to the technique 
discussed in Sect. 12.3.21 and ensures that signal variations inside a low-resolution pixel are accounted for. The high-resolution noise- 
weighted map consistent with th e above formalism is constructed as the first step of the map-making carried out by the MADmap 
program (iCantalupo et al.ll2009l) . The matrix eigenvalue decomposition is done using a ScaLAPACK0 interface that allows efficient 
parallel eigenvalue decomp osition of large matrices using a divide and conquer algorithm. 

The ROMA code (iNatoli et al.ll200lt Ide Gasperis et al.ll2005l) is an implementation of the optimal GLS iterative map ma king 
specifically designed for Planck, but also successfully used on suborbital experiments such as BOOMERanG (Masi et al . 20061). To 
estimate noise covariance as in Eq. (fTTt we start by calculating row i of its inverse, N~ l ,by computing its action on the unit vector 
along axis i: 




(30) 




(31) 




(32) 



(N~ l )t = (A T N~ l A) • e, . 



(33) 



2 http : //crd . lbl . gov/~borrill/cmb /madcap/ 

3 http : //www . netlib . org/scalapack/scalapac 



c_home . html 
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This calculation is implemented as follows: i) projecting the unit vector into a TOD by applying A on ef, ii) noise-filtering the 
TOD in Fourier space; iii) projecting the TOD into a map by applying A T . By computing each column independently we reduce 
memory usage because we allocate memory only for one map instead of allocating memory for the full matrix. The computational 
cost of the full calculation is dominated by FFTs that are repeated as many times as the number of columns, hence the scaling can 
be expressed as: 

Nflops ~ O (^samples " log 2 Ofilter) " ^pix) - (34) 

Once the inverse noise matrix is assembled it is inverted in a fashion similar to the one described above. We note that the resulting 
covariance matrix has to be symmetric. Though this is not a priori ensured by the algorithm it is the case within expected numerical 
errors. To ameliorate any effect we symmetrize the result by averaging the matrix with its transpose. 

Alternately, we can compute the NCM column-by-column with help of multiple map-making-like operations, i.e., 

yt = (A T N~ 1 A)~ 1 - eu (35) 

where ei is a unit vector as defined above and j; stands for a column of N. We rewrite the above equation as 

(A T N- 1 A)-y i = e i , (36) 

and solve it using the standard PCG map-making solver. We note that in such a case there is no need to store a full inverse noise 
covariance matrix in a memory of a computer at any single time as the operations on the left hand side can be performed from right 
to left. As a result this approach can be applied also for high-resolution cases for which the direct method described above would 
not be any more feasible. 

We note that unlike in the previous approaches based on the direct matrix inversion in the latter case there is no special car e 
taken of potential singularities. Though the presence of those does not hamper the PCG procedure (e.g.. ICantalupo et all [2009), 
nonetheless care must be taken while interpreting its results. 



3.2. Destriped map covariance 

In the destriping approach to map-making, we model all noise correlations by baseline offsets. Thus we write Eq. © as 

d = As + Bx + fi u , (37) 

where n u is a vector of uncorrelated white noise samples. Accordingly, we must replace the time-domain noise covariance matrix, 
N, by a diagonal matrix, N u . All noise correlation is then included in the prior baseline offset covariance matrix, V. 

If we now apply the destriping approximation to Eqs. (ITTI) we find for the pixel-pixel residual noise covariance matrix: 

M~ l = A T N~ l A - A T N~ l B + B T N~ l B)~ l B T N~ l A . (38) 

The first term on the rhs is the binned white noise contribution (a diagonal or block-diagonal matrix for temperature-only and 
polarized cases respectively) and the second term describe s the pixel-pixel correla tions due to errors in solving for the baselines, 
i.e., the difference between the solved and actual baselines (iKurki-Suonio et al1l2009|) . 

When making a map using destriping, one can use high resolution to solve for baselines and still bin the map at low resolution. 
Since this is equivalent to producing first a high-resolution map and then downgrading through inverse noise weighting, we will 
always assume the same pixel size for both of these steps. 



3.2.1. Conventional destriping 

Springtide (lAshdown et al.ll2007b|) is an implementation of the conventional destriping approach which solves for one baseline per 
pointing period. Since the baselines are so long, it allows for a number of optimizations in the handling of the data. During one 
pointing period, the same narrow strip of sky is observed many times, so the time-ordered data are compressed into rings before 
doing the destriping. Another effect of the long baselines is that the prior covariance matrix of the baselines, P, is strongly diagonal- 
dominant, so to a very good approximation can be assumed to be diagonal. As a consequence, the matrix (p~ l + B T N~ l B^j that 
appears in the expression for the inverse map covariance matrix (l38t is also diagonal. Thus the number of operations taken to 
compute (l38t is 

Nflops ~ 0(ftbase(ftpix/base) 2 )- (39) 

The number of baselines is small compared to the generalized destriping approach, so another method of computing the noise 
covariance matrix of the map becomes feasible. It is possible to compute the inverse posterior covariance matrix of the baseline 
offsets explicitly, to invert it and use it to compute the map covariance matrix. This method has the advantage that the resolution 
at which the destriping is performed is not constrained to be the same as the resolution of the final map covariance matrix. The 
destriping can instead be done at the natural resolution of the data to avoid subpixel striping effects. The inverse of the posterior 
baseline error covariance matrix can be calculated using Eq. (TT3T) and the corresponding map covariance matrix is given by Eq. (fT2l) . 
However, the pointing matrices, A, need not be for the same resolution in both steps. 

The structure of the inverse posterior baseline covariance matrix, Eq. (fT3K depends on the scanning strategy, but it is in general 
a dense matrix. Inverting the matrix and using it to compute Eq. (Tl2t involves dense matrix operations, so this method is computa- 
tionally more demanding than the other approach described above. However, the posterior baseline covariance matrix needs only to 
be computed once and stored, and then can be used many times to compute the map covariance matrix for any desired resolution. It 
is also possible to use Eq. (Tl2t to compute the residual noise covariance matrix for a subset of the pixels in the map. 
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Table 1. Pixel side to baseline length at 1 rpm spin rate 



side 


-y/pixel area 




^base 


4 


14°658 


192 


2.443s 


8 


7°329 


768 


1.222s 


16 


3!665 


3,072 


0.611s 


32 


1!832 


12,288 


0.305s 


64 


(£916 


49,152 


0.153s 


128 


(£458 


196,608 


0.076s 



3.2.2. Generalized destriping 

Madam (Keih anen et al J 120051: iKeihanene t al. 2009) is an implementation of the generalized destriping principle. It is flexible in 
the choice of baseline length and makes use of prior information of baseline covariance (P is not approximated to be diagonal). 
Even for a generalised destriper, all the matrices in Eq. (l38t are extremely sparse. Most of the multiplications only require 

operations proportional to the number of pixels, baselines or samples. The matrices P~ l and (P~ l + B T N~ l B > j are approximately 
circulant, band-diagonal matrices whose width is determined by the noise spectrum. For all cases studied in this paper, the latter 
matrix is limited to the order of 10 3 non-negligible elements per row. We call this width the baseline correlation length, n corr , and 
includes the white-noise contribution as well. n corr corresponds to a few hours of samples and is inversely proportional to baseline 
length, ftbi- 

We evaluate the prior baseline offset matrix from the power spectral density of the correlat ed noise, P(f), by Four ier transforming 
the baseline PSD, P x (f), into the autocorrelation function. The baseline PSD is evaluated as (Keihane n et alJl2009|) : 

1 °° sin^ 7ix 

PAf) = ~ — X P(f + ^/WWW + m), where g(x) = (40) 

'base m tl W 

The sum converges after including only a few m around the origin. For stationary noise, any row of P~ l can then be evaluated as a 
cyclic permutation of T~ l [1/P x (f)]. 

We evaluate (f38l) after computing (|40l) by approximating the inner matrix, (p- 1 +b t n- 1 b)~ , as circulant. This allows us to 
invert the ^base X ^base band diagonal matrix by two short Fourier transforms. It turns out that the matrix multiplications are most 
conveniently performed by first evaluating the sparse A T N~ l B matrix and then operating with it on the inner matrix from both sides. 

In effect the inverse covariance matrix gains a contribution from all quadruplets (jc,-, Xj, p, q), where baselines Xi and Xj are within 
baseline correlation length and hit pixels p and q. The number of operations required to complete the estimate is then proportional 
to 

^base " ^corr ' Opix/base) 2 ), (41) 

where ^base and n corr , the number of baselines per survey and correlation length respectively, are inversely proportional to the length 
of a baseline. In contrast, ft p i x /base, pixels per baseline, is proportional to baseline length. For short baselines and low-resolution maps, 
the magnitude of ft p i x /base is close to unity. Table [T] lists low-resolution A^e parameters and the baseline lengths that correspond to 
average pixel sizes. It can be used to estimate ft p i x /base and shows that, for example, 1.25 s baseline offset at A^de = 32 resolution 
covers approximately 4 pixels. The success of this approach is to replace n?, in the computation complexity by n 1 



pix/base' 



3.3. Smoothed covariance matrices 



In order to apply the smoothing operator to the low-resolution noise covariance matrices, N, (Eq. (l29t), we assume that low- 
resolution maps, and thus also their covariance matrices, are expected to cover the entire sky. We can also use the eigen- 
decomposition of the noise covariance matrices as it is available from the matrix inversion procedure described earlier. Consequently, 
we perform the smoothing to the eigenvectors of the noise covariance matrix, Eq. ([29b . 

Each eigenvector is a 3Af p i X map itself (I, Q, U map) and it has therefore an expansion in the spherical harmonic domain with a 
set of coefficients a t 

at = Ye t . (42) 

Here Y is the matrix that performs the transformation. It is made of spin-0 and spin-2 spherical harmonic functions. The coefficients 
shall now be smoothed with the same window function that we used on the high-resolution map: 

a t = Wa t . (43) 
Next we turn the smoothed coefficients back to a map 

% = Y-% (44) 
and compose the smoothed matrix 

N = Y J ^e i -eJ. (45) 
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We note that because N is symmetric, its eigenvalues are real and its eigenvectors make an orthogonal system. However, Ai and e\ are 
not in general the eigenvalues and eigenvectors of the smoothed matrix N. This is particularly important whenever the unsmoothed 
noise covariance matrix, N, is singular and therefore its calculation needs to be regularized, as described in Sect. 13. 1[ In such cases, 
special care may need to be taken to account for effects of such singularities on the smoothed covariance. As we point out in the 
next Section, the unsmoothed covariance is indeed commonly expected to be singular or nearly so and therefore a general procedure 
of treating singular cases is needed. We will discuss a proper way of dealing with such an issue in Sect. 16.21 

In addition, the smoothing procedure on its own will often lead to singular eigenvectors of the smoothed covariance matrix, 
with the eigenvalues corresponding to those close to zero. Though at a first glance such eigenvectors may look like being strongly 
constrained by the data, their actual value in the analysis is negligible as the sky signal in those modes is also smoothed. The 
nearly vanishing variance of those modes will often spuriously exaggerate the smoothing and map-making artifacts likely present 
whenever any inverse noise weighting needs to be applied. To avoid such problems hereafter we compute the inverse of the smoothed 
covariance via its eigenvalue decomposition and set the eigenvalues of all the nearly singular modes to zero. The criterion for 
selecting the nearly singular modes will in general depend on the case at hand. 

In some cases the eigenvalue decomposition of the smoothed matrix may not be readily available or its computation not desirable. 
We can then regularize the inversion of N by adding some low level of the pixel-independent uncorrelated noise. For consistency, 
a random realization of such noise should also be added to the corresponding maps. We note that both approaches are effectively 
equivalent and that the choice of the singularity threshold needed to select the singular eigenmodes corresponds roughly to the 
choice of the noise level to be added. We will commonly use the latter approach in some of the power spectrum tests discussed later. 



3.4. Singularities 

As we have pointed out in Sect. 13.11 the inversion of the inverse noise covariance matrix, N~ l , Eqs. (f20l) & (ITTT) . often needs to be 
regularized due to the presence of singular or numerically singular modes. In this Section we discuss the origin of such modes. 
We first note that in all cases considered here the inverse covariance can be expressed as 

AT 1 = A 1 MT l A, (46) 

where AT 1 is defined to be, 

{ TV" 1 , for the optimal maps; 

M = | TV; 1 - N~ l B + B T N- l B)~ l B t N~\ for the destriped maps. (47) 

We assume that the pointing matrix, A, has full column rank, and thus Ax = => x = 0. This is equivalent to an assumption that the 
sky signal can indeed be estimated from a given data set. Though this may not be always the case, in particular for the polarization 
sensitive experiments, it can usually be achieved if some of the ill-constrained pixels are removed from consideration. Given this 
assumption, the problem of the singular modes of N~ l becomes that of the matrix, Al" 1 , defined above. 

Let us consider the optimal map case first. The matrix At -1 is equal to the inverse of the time-domain noise covariance, N~ l . 
The latter, Sect. 12.21 is a block Toeplitz matrix with each block defined by an inverse of the noise power spectrum, Eq. (|2T1) . Each 
of those blocks describes the noise properties of one of the stationary data segments assumed in the simulations. For each block the 
eigenmodes corresponding to the lowest frequencies as permitted by the length of the segment have eigenvalues vastly smaller than 
the high frequency modes. These modes can therefore lead to near singularities. This is specifically true for zero-frequency modes 
corresponding to an offset of each of the stationary data segments. The (near) null space of the full matrix will be therefore spanned 
by all such vectors corresponding to each of the segments. 

Due to projection effects, not all of those modes result in singular modes of the final pixel domain noise covariance. However, if 
for a mode, t, from the null space of TV -1 , there exists a pixel domain vector, jc, such as t = A x, the inverse noise covariance will 
be singular with eigenvector equal to x. 

In the studied case, the scanning strategy is such that the sky areas observed in each of the stationary periods overlap, which 
efficiently removes most of the potential degenerate vectors. In fact only a single pixel-domain IQU vector of which the / part is 
one and all the others zero, called hereafter a global offset, can potentially be singular. We will indeed confirm these expectations 
via numerical results later 0. 

In the case of the maps produced with the destriper codes in the absence of any priors, i.e., !P _1 = the M~ l matrix has as 
many singular vectors as the baseline offsets defined by the columns of the offset 'pointing' matrix, B. However, as long as all of 
the offsets cross on the sky the only singular vector of the pixel-domain covariance will again correspond to the global offset vector 
as in the optimal map case. We note however that unlike in that case, this time this vector is exactly singular. If a prior is employed, 
as is the case in both the classical and generalized implementations of the destriper technique discussed here, the columns of the 
matrix B are no longer singular vectors of the matrix At -1 , nor is the global offset vector a singular vector of N. Nevertheless, at 
least for some common choices of the prior the global offset vector remains nearly singular. 



4 We note that in the argument presented here the global offset mode is only nearly singular. This is due to our assumed noise power spectrum, 
which is finite at zero frequency. In a more realistic case the offsets of the stationary segments will be however unknown, corresponding to an 
infinite amplitude at zero frequency. The noise covariance in such cases should be therefore considered to be strictly singular with the global 
offset being the singular eigenvector. In such cases the noise weighting on the right term of the map making equation, Eq. [T9j will force the offset 
of each of the stationary data segments to be strictly zero. This may not be a sufficiently good approximation in particular for short stationary 
time segments. This could, however, be alleviated by introducing the segment offsets as extra degrees of freedom contained in the vector, x, (e.g. 
IStompor et alll2002h . 
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4. Numerical tests and comparison metric 

This paper has two main goals. On the one hand we propose and compare various methods devised to produce the low-resolution 
maps, searching for the map-making method which leads to low-resolution maps virtually free of artifacts such as those due to sub- 
pixel power aliasing. In parallel we develop the tools to estimate residual noise covariance for such maps, which properly describe 
the error of the estimated maps due to residual noise. The numerical results presented in the subsequent sections of the paper aim 
therefore at comparing and validating the algorithms which we have presented earlier. The discussed comparisons involve standard 
statistical tests, such Kolmogorov-Smirnov,^ 2 , etc. and ones which are specifically devised in the light of the anticipated future 
applications of the maps and their covariances. We describe such tests in this Section. 



4.1. Quadratic maximum likelihood power spectrum estimation 

One of the main applications for the low-resolution noise covariance matrices we envisage is to the estimation of power spectra, Q. 
This is often separated into the estimation of large and small angular scales, usually associated with high and low signal to noise 
regimes. The methods discussed here are relevant for large angular scales. The successful estimation of the underlying true power 
spectrum of the sky signal sets demanding requirements for the quality of the maps produced for such a purpose as well as the 
consistency of the estimated noise covariance and the actual noise contained in the map. For this reason we will use hereafter power 
spectrum estimation as one of the metrics with which to evaluate the quality of the proposed algorithms. 

W e will use the Quadratic Maximum Lik elihood (QML) method for the powe r spectrum estimation as introduced in ( Tegmark 
1997) and later extended to polarization in (iTegmark & de Qliveira-Costall2QQll) . Given a map in temperature and polarization 
m = (T, Q, 17), the QML provides estimates of the power spectra, that read, 

$ = Z (F_1) ^ r W E x>™ - tr(A^)] • (48) 

Here, Cf is an estimated power spectrum, X = TT, EE, TE, BB, TB, or EB, and F xx , is the Fisher matrix defined as 
1 



C -,0C C _ X 8C 



acf acf 

The E matrix is given by 



(49) 



E x =\ C - l ^-C-\ (50) 
x 2 dC* 

where C = 5(Q) + N is the covariance matrix (signal plus noise contribution) of the map, m. Here Q is a fiducial power spectrum 
needed for the calculation of the signal part of the covariance. In this paper we will take it to be given by the true power spectrum 
as used to produce the simulated skies. Though this would be an unfair assumption, while testing the performance of this power 
spectrum estimation technique, it is justified in our case, where the fact that it leads to the minimal estimation uncertainties increases 
the power of our test. (Indeed the QML estimator is in fact al so known to be equivalent to a single iteration of a quasi-Newton- 

Raphson procedure to search for the true likelihood maximum dBond et al.ll 19981).) 

M ore details about t he QM L method can be found elsewhere (e.g. lTegmarkl 19971: ITegmark & de Oliveira-Costal200ll; lEfstathiou 

2006). Gruppu so et al.l (120091) describes the specific implementation of the method, nicknamed Bolpol, as used in this work and 
discusses its performance in the application to the WMAP 5 year data. 

Hereafter, we neglect any systematic effects of either instrumental and/or astrophysical origins. Nevertheless we note that if 
the final CMB map is obtained via some linear cleaning procedure involving maps computed either for different detectors and/or 
frequency channels, the results of this paper will be still relevant and the noise covariance of the 'cleaned' CMB map can be 
computed via a linear combination of the single map covariances calculated in turn with help of the procedure discussed here. 



4.2. Noise bias 

In the map-making methods considered in this paper residual noise in the maps is independent of the sky properties and completely 
defined by the time-domain noise properties and the scanning strategy. The noise present in the maps contributes to the power 
spectrum estimates of the map signal. We will therefore refer to this contribution as the noise bias and use it to quantify the noise 
level expected in the maps of our different methods in a manner more succinct and manageable than the full noise covariance. 
The noise bias is defined to be the expectation value for angular power spectrum in the noise-only case 

N f Y = (Cf' noise >, (51) 

where X and Y stand for T, E and B. In the general case of a quadratic estimator, such as QML, Eq. (|48lh 
(iTegmark & de Oliveira-Costall200ll) . the power spectrum estimates are given as a quadratic form of the input map, m, 

Cf Y = m T Qf Y m = tr Q$ Y mm T . (52) 

By taking m to be noise-only, the noise bias can be expressed as 

Nf Y = (Cf Y > noisG ) = trQX Y N. (53) 
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Table 2. Frequently used symbols in this paper. 



Symbol 


Definition 


C 


pixel-pixel covariance matrix 


N 


noise covariance, map domain 


N 


noise covariance, time domain 


AT 


correlated (1//) noise covariance, time domain 


M~ l 


white noise covariance, map domain 


K 


white noise covariance, time domain 


M p 


3x3 observation matrix 


V 


prior baseline offset covariance 


A 


detector pointing matrix 


B 


offset-to-TOD matrix 


X 


baseline offset vector 


m 


3A^ pix Stokes I,Q,U map vector 


s 


map estimate 


n 


noise vector 


n' 


correlated noise vector 


n u 


white noise vector 


d 


TOD vector 


s 


sky map 


e t 


Cartesian unit vector along the i:th axis 



Given the eigenvalue decomposition of the noise covariance matrix, N = 2i we can evaluate the noise bias as 

Nf^^ejQfe, (54) 

i 

In the following we will validate the noise covariance matrices, estimated for the recovered sky maps, with the help of their respective 
noise biases, which we will compare to results of Monte Carlo simulations. In this context it is particularly useful to consider a 
pseudo-Q estimator that assumes uniform pixel weights and full sky. For this estimator, the operator Q* Y is 

Q * ¥= 2eTiW Y *> (55) 

where has 2^+1 rows that are maps of the appropriate spherical harmonics. It maps a map vector into a vector of spherical 
harmonic expansion coefficients {of } where m = -£ . . . £. 

The procedure we implement here involves two steps. First for every estimated noise covariance we compute analytically 
the noise bias using Eqs. (l54l) and (l55t . Second, we compute the bias using Monte Carlo realizations of the noise-only maps pro- 
duced using the cor responding map-making procedure. The multiplications Yfui are conveniently implemented using the HEALPix 
(IGorski et al.l l2005) Fortran 90 subroutine map2alm. We note that as we consider hereafter only full-sky cases, the noise biases we 
compute as described above would be equal to those expected in the Maximu m Likelihood estim ates of the sky power spectrum, 
were it not for the imperfection of the sky quadrature due to pixelization effects (IGorski et al.ll2005l) . 

5. Simulation 

5.1. Scanning strategy 

In this study the Planck satellite orbits around the second Lagrangian point (L2) of the Earth-Sun system (Pu pae & Tauberl feOCM). 
The spin axis lies near the ecliptic plane, precessing around the anti-Sun direction once every six months with an amplitude of 7° 5. 
The telescope line-of- sight forms an 85° angle with the spin axis. In addition to these modes, we in clude a nutation of the spin axis 
and slight variations to the 1 rpm spin rate. Details of the scanning simulation can be found from lAshdown et afl (|2009|) where it 
was used in a map-making study. 

5.2. Planck detectors 

In this Paper we study residual noise in the Planck 70 GHz frequency maps. Planck has twelve detectors at 70 GHz. In the focal 
plane they are located behind six horn antennas, a pair of detectors ("Side" and "Main" detector© sharing a horn. A pair of 
detectors measures two orthogonal linear polarizations. The horns are split in two groups (three horns in a group). The Side and 
Main polarization sensitive axes of a group are nearly aligned and the polarization directions of the second group differ from the 
first group nominally by 45°. Two horns from the different groups make a polarization pair that follows the same scan path in the 
sky (three pairs in total with slightly different scan paths). As a minimum the observations of a polarization pair are required to 
build a polarization map. Due to implementation restrictions the Side and Main polarization axes are not fully orthogonal and the 



5 Side and Main refer to two detector branches downstream from the orthomode transducer that separates the two orthogonal linear polarizations. 
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polarization direction differences between the groups are not exactly 45°, but the deviations from these nominal values are small 
(< 0?2). The Side polarization axes of the two groups differ by +22?5 and -22?5 from the scan direction. 

The beams of the detectors were assumed circularly symmetric with a 14' FWHM (full width half maximum) beam width. 
The beams do not impact the residual noise maps or covariance matrices. None of the map-making methods studied here make an 
attempt to correct for beam effects in the maps. 

5.3. Time ordered data 

We computed the NCM's of our three map-making methods using our noise model spectrum and the one year pointing data of 
twelve 70 GHz detectors. We produced NCM's for both N^q = 8 and N^q = 32 pixel sizes. We wanted to compare these NCM's 
to the noise maps made by the same map-making methods. For that purpose we simulated 50 noise-only timelines and made maps 
from them. Ou r correlated noise streams were simulated in six day chunks by inverse Fourier transforming realizations of the 
noise spectrum (iNatoli et al .1120021) . We assumed an independence between the chunks and between the detectors. Fig. Q] contains a 
comparison between the power spectra of the generated noise streams and the model spectra. 

Twenty-five of the surveys featured a relatively high 1 // contribution having the knee frequency, /knee, set to 50 mHz. The other 
half was simulated to have a more favorable /knee = 10 mHz. It should be noted that these frequencies have been chosen above and 
below the satellite spin frequency, 1 rpm « 17 mHz. The slope of the 1// noise power spectrum was a = -1.7. The correlation 
timescale of the 1 // noise was restricted to about one day. This made our noise spectrum flat at low frequencies (below a minimum 
frequency / min = 1.15 x 10" 5 Hz * 1/24 h). As we described earlier, it is the minimum frequency that determines the correlation 
length of the noise filter in the optimal map-making. In the noise covariance matrix of the generalized destriping, the baseline 
correlation length is, however, determined by the knee frequency. 

We used a uniform white noise NET of 204 iaK Vs for all detectorfj. We chose this NET because we wanted to produce noise 
maps and covariance matrices whose noise levels are compatible with another CTP study dJaffe et al .1 12009). In all map-making and 
NCM computations we assumed a perfect knowledge of the detector noise spectrum. 

The noise timelines were processed directly into both low-resolution (7V S id e = 32) and high-resolution (7V S id e = 1024) HEALPix 
maps using the discussed map-making codes. In the A^de = 1024 temperature and polarization maps the mean standard deviations 
of white noise per map pixel were 44 and 63 yuK (Rayleigh- Jeans yuK). For A^e = 32 maps the corresponding values were 1.4 and 
2.0//K. 

The high-resolution maps were in turn downgraded to the low target resolution using the schemes detailed in Sect. 12.31 
For the signal error studies described in Sect. 16.31 we scanned simulated foreground maps into signal-only timelines. These 
we processed into low-resolution (7V S id e = 8) maps using the same methods as for the Af s id e = 32 (both directly and through high 
resolution). We then extracted the signal error part by subtracting a binned map from the destriped map. The foreground signal 
errors were summed with a CMB map to provide a worst case scenario of signal striping in otherwise perfectly separated CMB 
map. 

5.4. Input maps 

To study bandwidth limitation with respect to downgrading we simulated 117 high-resolution Af s id e = 1024 CMB skies corresponding 
to the same theoretical spectrum, Q. These maps were smoothed and downgraded to /V S ide=8 using three different Gaussian beams 
of widths 5°, 10°, 20° and three apodized step functions with the choices of (£i,£ 2 ) being (20, 24), (16, 24) and (16, 20). A seventh 
set of downgraded maps was produced by noise weighting according to the scanning strategy. To comply with this last case, the 
smoothing windows include also the /V S ide=8 pixel window function from the HEALPix package. 

For the signal error exercise we used the Planck sky model, PSJVfl version 1.6.3, to simulate the full microwave sky at 70 GHz. 
For diffuse galactic emissions we included thermal and spinning dust, free-free and synchrotron emissions. We then added a 
Sunyaev-Zeldovich map and finally completed the sky with radio and infrared point sources. The combined Af s id e =2048 map was 
smoothed with a symmetric Gaussian beam and scanned into a timeline according to the scanning strategy. 

For final validation the noise covariance matrices were tested in power spectrum estimation. Each noise map was added to a 
random CMB map drawn from the theoretical distribution defined by a fixed theoretical CMB spectrum. The theoretical spectrum 
is the WMAP first-year best fit spectrum and has zero BB mode. 

All maps in this work are presented in the ecliptic coordinate system. This choice is useful for Planck analysis since the scanning 
circles and many map-making artifacts form circles that connect the polar regions of the map. In this coordinate system the galaxy 
is not positioned in the ecliptic plane but forms a vertical horse shoe shape around the center of the map (see Fig.[T6l). 

6. Results 

In this Section we first focus on the noise covariance matrices computed using different map-making techniques. We discuss and 
compare the overall noise patterns implied by such matrices and test the quality of the destriper approximation as applied to the noise 
covariance predictions. In the second part of the Section we discuss the low-resolution maps and evaluate their quality in the light 
of their future potential applications, such as those to the large angular scale power spectrum estimation. Due to the computational 
resourc e res trictions the low-resolution results presented here a re obtained either with the HEALPix N^q = 8 — in particular tests 
in Sect. l6.3l and power spectrum estimation tests in Sect. 16.2.31 — or Af s id e = 32 — in most of the other Sections. 

6 <t = NET y /sample, where <x and /sample were defined in Eq. (12TT) . The 70 GHz detectors had /ample = 76.8 Hz. 

7 see http : //www . ape . univ-paris7 . fr/APC_CS/Recherche/Adamis/PSM/psky-en . php 
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Figure 2. Eigenspectra of the inverse covariance matrices N~ l . 
MADping, ROMA and Madam results for /k^ =50 mHz over- 
lap completely. Springtide results are for lOmHz. 



6.1. Noise covariance matrices 

First we discuss the noise covariance matrices computed for the low resolution maps of the direct method. As explained in Sect. I2.3I 
we compute from these matrices the noise co variances of the other downgrading techniques. In the following we consider noise 
covariance calculated using 4 different ways. In the first way we compute the noise covariance using the optimal algorithm. For 
this purpose we have developed two codes MADping or ROMA, which are described in Sect. I3.ll However, as they are just two 
different implementations of the same algorithm, we derive most of the results presented in the following and involving the optimal 
covariance using MADping. We note that whenever results from the both codes are available they have turned out to be virtually 
identical within the numerical precision expected from this kind of calculations. The optimal noise covariance matrices are expected 
to provide an accurate description of the noise level found in the actual optimal maps. We will test this expectation in the following 
and use the optimal results as a reference with which to compare the destriping results. 

The three remaining computations of the noise covariance are based on the destriping approach and correspond to different 
assumptions about the offset prior as well as baseline length. We consider the following specific cases: a classical destriper calcula- 
tion with a baseline of 3600 s (Springtide) and two generalized destriper computations with a baseline of 1.25,s and 60 s (Madam). 
For each of these cases we will compare the covariance matrices with each other, with the optimal covariances and then test their 
consistency with the noise found in the simulated maps. We note again that this last property is not any more ensured given the 
approximate character of the destriper approach. 

Fig. [2] shows the eigenvalue spectra of some inverse NCMs. We note that all matrices possess a positive semi-definite eigen- 
spectrum as is required for any covariance matrix, yet at the same time they all have one nearly ill-conditioned eigenmod^l , which 
renders the condition number, i.e., the ratio of the largest and smallest eigenvalue, very large. This is in agreement with our expec- 
tations as described in Sect. I3.4I Indeed the peculiar eigenmodes corresponding to the smallest eigenvalues of the inverse matrices 
are also found to be non-zero and constant for the / part of the vector and nearly zero for its polarized components, and thus close to 
the global offset vector discussed in Sect. 13. 41 The small deviations, on the order of 10~ 3 , exist, as expected, as none of the peculiar 
modes is in fact truly singular. We note that the MADping and ROMA results, both computed in this test, are seen in the figure to be 
indistinguishable. They also coincide very closely with the Madam results computed for the same /knee = 50 mHz. The Springtide 
results, computed with /knee = 10 mHz are close to, though not identical with, the Madam results for the very same value of /k nee 
when a longer (60 s) baseline is used for Madam. 

Fig. [3] depicts the estimated Stokes I, Q and U pixel variances as well as the covariance between them. These quantities are 
dominated by the white noise contribution and all methods describe white noise in the same manner. The top right-most panel 
shows the reciprocal condition number (1/condition number) of the 3 x 3 blocks of the matrix A T N~ 1 A for each of the sky pixels. 
These numbers define our ability to disentangle the three Stokes parameter for each of the pixels. Whenever they are equal to 1 /2 
the parameters can be not only determined but their uncertainties will not be correlated. If the reciprocal condition number for a 
selected pixel approaches 0, the Stokes parameters can not be constrained. In the cases considered here, the Stokes parameters can 
clearly be determined for all the pixels. 

Fig. [3] shows a strong asymmetry between the IQ and IU blocks. The fact that the polarization axes of the two detectors of a 
horn are not fully perpendicular makes the I noise of a pixel correlate with the Q and U noises of the same pixel. We can imagine an 
instrument basis, where one group of three horns measures the Q of this basis and the other group measures the U of the same basis. 
Because Side and Main polarization axis deviations from the orthogonality have similar magnitudes in the two groups, we expect 
the diagonals of IQ and IU blocks to be similar (symmetric in IQ and IU) in the instrument basis. In the map we use a different 



8 A double precision (64 bit) matrix is numerically ill-conditioned when the condition number exceeds approximately 10 12 (Press 1992). For 
our matrices the situation is not as dire but the first eigenmode still deserves special attention. 
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polarization basi^J Building the noise weighted map (in the map-making) rotates the Q and U from the instrument basis to the map 
basis. The IQ and IU blocks become asymmetric in this rotation. 




Figure 3. Top: MADping pixel variances for temperature and polarization and the reciprocal condition number of the pixel observa- 
tion matrices. Bottom: Correlation coefficient xlO 3 between I-Q, I-U and Q-U pixels. This part of the noise covariance is dominated 
by white noise which is modelled equivalently in all three paradigms. Hence, MADmap, ROMA, Madam and Springtide results are 
nearly identical. Maps are rotated into galactic coordinates to show the structure near the ecliptic poles. 

Figs. 0]-[5] show plots of a single column of the noise covariance matrix. The column corresponds to reference pixel number 
0. In the HEALPix nested pixelization scheme for A^e = 32 resolution, pixel is located at the equator. In the plots each pixel 

has the value (m p m q ) normalized by ^(m 2 )(ra 2 ). Thus the pixel values of the plots represent correlation coefficients. Due to this 

normalization, the reference pixel automatically gets unit value and is later set to zero in order to bring out smaller features of the 
other elements of the columns. 

The sky is scanned from one ecliptic pole to the other. The NCM column maps are characterized by bands of correlation along 
the scanning rings. Pixels near the equator, such as the reference pixel 0, are only observed during a few-hour window as the 
satellite scanning ring is rotated over the course of the survey. The two crossing bands of higher correlation correspond to two 
pointing periods half a year a part that observe the reference pixel 0. 

For both generalized destriping and optimal map-making, there is a visible gradient in the correlation along the scanning ring. 
Pixels that are observed immediately before or after the reference pixel have the strongest correlation. The conventional destriping 
with its hour long baselines assumes constant correlated noise over the scanning ring and does not, therefore, show this feature. 

Figs. @] and [5] show that the strongest cross-correlations (~ 1%) exist between Q and U noise maps. Side and Main polarization 
sensitive directions differ slightly from 90° and, as a result, small (~ 10" 4 %) IQ and IU correlations remain. 

Figs. [SflT] show plots of the NCM columns of another reference pixel (number 2047). This pixel is located at the northern ecliptic 
pole region and it exhibits a very different correlation pattern compared to the previous case. Since the pole is visited frequently 
through the course of the survey, it becomes correlated with the rest of the map as a whole. Correlation amplitude is increased by a 
factor of 3 from the equator pixel case (the increase can be seen from the color bar ranges) and there is now a distinct asymmetry 
between northern and southern hemispheres. As one expects, the asymmetry only appears in optimal and generalized destriping 
estimates. 

6.2. Noise covariance validation 

In this Section we report the results of three different validation tests. We performed a^ 2 test, compared the noise biases computed 
from the matrices and corresponding Monte Carlo maps and finally used the matrices and Monte Carlo maps as inputs to angular 
power spectrum estimation. 

6.2.1. By X 2 

Residual noise expected in the recovered maps is Gaussian due to the linear character of all the map-making methods considered 
here. Thus the noise is completely described by its covariance matrix. More specifically, in the absence of any singular modes of 



9 In a HEALPix map the Stokes parameters Q and U at a point in the sky are defined in a (x,y,z) reference coord inate, where the x-a xis is along 
the meridian and points to south, the y-axis is along the latitude and points to east, and the z-axis points to the sky (Gorski et al. 2005). 
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Figure 4. A single column of the MADping covariance matrix corresponding to a pixel at the ecliptic equator. For both ROMA and 
Madam 1.25 s baseline counterparts all visible characteristics remain unchanged. We plot the value of the correlation coefficient, R, 
multiplied by 10 3 . In order to enhance the features, we have halved the range of the color scale. 
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Figure 5. A single column of the Springtide covariance matrix corresponding to a pixel near the ecliptic equator. For description of 
the normalization, see text. 
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Figure 6. A single column of the MADping covariance matrix corresponding to a pixel at the ecliptic pole. For both ROMA and 
Madam 1.25 s baseline counterparts all visible characteristics remain unchanged. We plot the value of the correlation coefficient, R, 




Figure 7. A single column of the Springtide covariance matrix corresponding to a pixel near the ecliptic north pole. For description 
of the normalization, see text. 
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the estimated residual covariance, N, the residual maps, m = s - s, are drawn from a multivariate Gaussian distribution defined by 
N. Therefore, the^ 2 statistic, defined as,;^ 2 = m T N~ l m, is drawn from a^ 2 distribution with 3N V [ X degrees of freedom (dof). 

If any singular mode is present we simply replace the matrix, iV -1 , in the definition of x 2 , by a matrix, iV' -1 , which is like N~ l 
in all respects but has the eigenvalue corresponding to the singular vector set to 0. We note that if the eigenvalue decomposition of 
the matrix iV -1 is not available or too costly to compute, we can achieve numerically the same effect by defining N' = N + rf vv\ 
where v is a sing ular vector we want to project out and rj 1 is a large positive number for which however inversion of N' is still stable 
fe.g. lBond et al]l200Qh . As we subtract one degree of freedom corresponding to the excluded, ill-conditioned eigenmode we expect 
that there are in total 3N P i X - 1 degrees of freedom left in our maps 

We can apply the analogous test to the smoothed noise covariance. The x 2 statistic is defined as before with the smoothed 
covariance matrix as well as residuals used now in place of the respective unsmoothed objects. As we commented on that in 
Sect. 13.31 the inverse of the smoothed covariance has to be appropriately regularized, to avoid the results of the test being biased 
by the artifacts potentially present at the scales smaller than the smoothing kernel and therefore not containing any cosmologically 
useful information. The effective number of degrees of freedom left in the data will coincide then with the number of eigenvalues 
which have not been set to zero in the regularization process. Alternately, if the preferred regularization approach involved adding 
some low level of the white noise, Sect. 13.31 the number of the degrees of freedom is equal to 3N V i X . 

In addition, one needs to take care of the singularity of the unsmoothed noise covariance. This has to be done explicitly if 
the smoothed version of the ill-conditioned eigenmode, v, does not belong to the null space of the inverse smoothed NCM, i.e., 
N~ l (Lv) ¥ 0. To do so, we employ the same approach as before, replacing the regularized inverse of the smoothed covariance 
matrix, AT 1 , by, 

AT 1 -> [ft + t] 2 Lv (Lv) T ] _1 AT 1 - (AT 1 Lv) [(Lv)' AT 1 (Lv)]" 1 (Ar 1 Lvf (56) 

where L is a smooth ing operator, Eq. ([29b , and the last expression follows from the Sherman-Morrison- Woodbury for- 
mula ( Woodbury 1950). This last operation additionally reduces the number of degrees of freedom by 1 (or whatever number 
of modes, v is to be projected out). 

The Kolmogorov-Smirnov test can be used to test whether a set of samples conforms to some theoretical distribution. The test 
estimates the probability of the maximal difference between the empirical distribution function, 

1 n 

F n (x) = -Ye(xi), (57) 
n 

of the observations jc,- (in our case the individual^ 2 ) and the theoretical cumulative distribution function. @(x) is the Heaviside step 
function. We note that in this work we take an advantage of the fact that we can simulate the residual noise directly. Though this is 
clearly not the case when real data are considered, the tests described here can be applied to a difference of two sky maps produced 
by dis joint sets of detecto rs operating at the same frequency and can therefore be a useful test of the real life data processing integrity 
(e.g. lStompor et al.ll2QQ2b . 

Figs. ISUTOl show the two cumulative distribution functions for 25 noise maps in the case of the direct method. Reported p- 
values are the probabilities of observing this level of disagreement even if the noise description was exact. Conventionally, the level 
p < 0.05 is considered to be enough to reject the null hypothesis that the distributions match. 

We then proceed to study the agreement between downgraded noise maps and the noise covariance matrices. As a test case, 
we use the Madam NCM for lOmHz knee frequency, 1.25 s baselines. We smoothed the covariance matrix using an apodized 
window function, setting the thresholds to 2Af S id e and 3Af S id e respectively. As expected, the smoothed matrix is extremely singular. 
We compute its inverse by including only the eigenvalues that are greater than 10" 2 times the largest eigenvalue, including 20, 882 
of the 36, 864 available modes. 

Fig. [TT] shows the empirical distribution functions of the^ 2 . Even though the matrix is computed for the direct method, the 
inverse noise weighted (INW) maps conform well to it. However, when we apply the smoothing kernel to the high-resolution maps, 
there is a clear disagreement. This stems from the fact that in this downgrading the high resolution pixels are not correctly (inverse 
noise variance) weighted when we compute the low-resolution map. If we first produce a low-resolution direct method or INW map 
and then smooth it, the agreement is much better. This is shown in the bottom row of Fig.fTTI 

In the case of the direct method maps our results show that only optimal (noise) maps and their respective noise covariance 
are mutually consistent in the light of the x 2 statistics. The good statistical agreement in this case does not depend on the time 
domain noise characteristics nor map resolution. This is expected given that the noise covariance estimator implemented in the 
optimal codes, Eq. (l20l) , is an exact expression describing the noise properties in the pixel domain, and that we have assumed 
perfect knowledge of the time domain noise. 

The level of consistency found in the destriping cases varies depending on the underlying time domain properties, i.e., /knee, and 
on the assumed baseline length. In the case of /knee = 50 mHz, we have not found a satisfactory agreement in any of the considered 
destriping cases. For the lower /kn ee = 10 mHz the results obtained with the generalized destriper, Madam, are satisfactory for the 
short, 1.25 s, baseline choice, and marginal for the long one, 60s. The latter result is consistent also with that obtained using the 
classical destriper, Springtide. 

In the case of direct low-resolution map-making, the discrepancy between the noise maps and the noise covariance matrix stems 
from the destriping approximations. Both destriping approaches assume that the correlated part of the noise is perfectly modelled 
by a set of baseline offsets. Correlated noise occurring at frequencies higher than what the baselines can model is not removed from 
the TOD and therefore is binned onto the map For short baselines or low knee frequencies the unmodelled noise manifests as a 
relatively small angular scale correlation and does not bias the power spectrum estimates at low I. 
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Figure 8. MADping empirical x 2 distribution 
function from the 25 residual noise maps com- 
pared with the theoretical cumulative probabil- 
ity density. The black stair line is the empiri- 
cal distribution function, the blue solid line is 
the theoretical x 2 distribution for 3A^ pix - 1 = 
36, 863 degrees of freedom. It is the same for 
all direct method maps in Figs. [8UT0l The red 
dashed line is the least squares fit of the^ 2 dis- 
tribution to the experimental distribution (dof 
being the fitting parameter). The horizontal axis 
is translated to the expected center of the distri- 
bution, (x 1 ) - dof = 36, 863, and scaled by the 



expected deviation, cr x 



/2dof. 



Figure 9. Madam empirical x 2 distribution 
functions from the 25 residual noise maps com- 
pared with the theoretical cumulative probabil- 
ity density. 



Figure 10. Springtide empirical^- 2 distribution 
functions from the 25 residual noise maps com- 
pared with the theoretical cumulative probabil- 
ity density. Left: direct method. Right: INW. 



We conclude that the noise covariance of the low-resolution maps produced by the destriping algorithm needs to be used with 
care. The flexibility of the generalized destripers permitting them to use different baseline lengths makes them in this context the 
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Figure 11. Madam empirical^ 2 distribution functions from 60 residual noise maps compared with the theoretical cumulative 
probability density. For a smoothed map, we count the degrees of freedom as the number of included eigenmodes in the inverse 
NCM. Top: Three sets of low-resolutions maps using just one of the low-resolution methods at a time. The high-resolution maps 
for INW and smoothing methods had A^ide = 1024. Bottom: Since it is suboptimal to compute a low-resolution spherical harmonic 
expansion from a noisy high-resolution map, we test how well the smoothing approach works in conjunction with the two pixel- 
based downgrading methods. 



preferred choice. We emphasize however that, if the accuracy of the noise description is the major concern, then only the optimal 
techniques are suitable. In the next Sections we will reconsider all these low-resolution map-making techniques in the context more 
specific to the large angular scale power spectrum estimation work, which is envisaged as the main application of the low-resolution 
maps and their covariances. 

6.2.2. By noise bias 

In this Section we describe the calculation of the average angular power spectrum of noise maps, i.e., the noise bias — see Sect. 14.21 
using a pseudo Q estimator for which both the NCM estimate and the Monte Carlo map averages are feasible to compute. Testing 
the noise covariance matrix by comparing estimated and measured noise biases can be viewed as complementary to the x 2 t ests 
described in the previous Section. It can certainly provide more information than the plain x 2 test, as instead of simple pass or fail 
indicator, the noise bias comparison will tell at which angular resolution the noise model agrees with the data. At the same time, 
the noise bias is less sensitive to the anisotropic features present in the residual noise. The noise bias test is clearly more directly 
relevant for power spectrum estimation. 

Figs. [T2l - fT4l compare noise bias averages from 25 noise realizations of the maps of the noise residuals to the analytical estimates, 
Eq. d54lh based on the estimated noise covariance matrix. Map spectra are computed using the HEALPix anaf ast utility and the 
estimated noise biases using the corresponding map2alm subroutine. We only show the autospectra, TT, EE and BB, as the scanning 
has decoupled the modes to large extent and only minimal coupling between the modes exist. 

The error band around the averages is the standard deviation of the individual Q values divided by V25. Each plot exhibits up 
to five curves: direct, noise weighted and harmonic smoothed noise biases, and two analytical estimates. 

The results derived for the case of the Madam runs with the short baseline, 1.25 s, and the high knee frequency, 50mHz, as 
shown in Fig. [T2] agree now very well with the MADping results. This is unlike in the chi 2 test discussed earlier, indicating that 
those were the anisotropic features responsible for the latter disagreement. The numerical calculations of the noise bias are in this 
case agree very well with the analytic predictions, Sect. 14.21 Both these facts validate the destriper approximation to the noise 
covariance in the light of this test, which is found to describe sufficiently precisely noise in the low-resolution map. This conclusion 
agrees with those derived using the^ 2 test earlier. 

The long baseline case is shown in Fig. [13] for = 50mHz and computed using the generalized destriper, and in Fig. [14] 
for the low value of the knee frequency, lOmHz, and based on Springtide results. We find that in the high knee frequency case 
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the prediction and numerical results differ rather dramatically, highlighting the failure of the destriper approximation in such case 
already seen with the x 2 t est - We note here that the failure seems to be affecting the largest angular scales as both the numerical 
and analytical results tend to converge at the highest £ end considered in this analysis. Similar results are found for the classical 
destriper maps and covariances. For the low /knee case the agreement is found to be marginal, with visible deviations seen generally 
at £ < 5 and are the most significant in the case of the BB mode spectrum. The results obtained using Madam in the analogous case 

are nearly indistinguishable. 

We note that our analyti c prediction are well i n line with WMAP findings ( Hins haw et al.ll2007l;|Page et al.ll2007l) and studies of 
the destriping framework ( Ef stat hioull20Q5[ 120061) . 

TT noise bias, 50mHz, 1.25s EE noise bias, SQxxiHz, 1.25s BB noise bias, SQixiHz, 1.25s 




Multipole, I Multipole, I Multipole, t 

Figure 12. Analytical and mean Monte Carlo noise biases from Madam runs at /knee = 50mHz using a short 1.25 s baseline. Grey 
band is the 1-cr region for the average, computed by dividing the sample variance by V25. 
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Figure 13. Analytical and mean Monte Carlo noise biases from Madam runs at /knee = 50 mHz using a long 60 s baseline. Grey 
band is the 1-cr region for the average, computed by dividing the sample variance by V25. Long baselines clearly fail to model the 
correlated noise. 
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Figure 14. Analytical and mean Monte Carlo noise biases from Springtide runs at /knee = 10 mHz. 

Averaged noise biases conform to the analytical estimates on almost all accounts and deviations are small compared to the 
absolute noise bias. As a warning example we include the results of modelling residual noise from the /knee = 50 mHz timelines 
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Table 3. Comparison of TT power (rms) in the stripe maps at A^de = 8 



Map-maker 


CMB a 
Direct 


Averaged 


Smoothed 


Foregrounds' 5 
Direct Averaged 


Smoothed 


Madam, 1.25s 
Madam, 60s 


5.44 //K 

2.07 yL/K 


7.29 nK 
3.31 nK 


5.20 nK 
1.98 nK 


9.46 fiK 4.67 nK 
3.81 yt/K 3.21 nK 


2.14nK 
2.23 nK 



a Binned CMB map rms is 37.34 /iK 

b Binned foreground map rms is 196.3 yt/K 



using 60 s baseline offsets. This leaves more noise in the maps but the corresponding analytical estimate is actually lower, since 
noise not modelled by the baseline offsets is neglected. 

At 10 mHz the Madam results for long 60 s baselines are equivalent with Springtide results: both succeed to estimate noise bias 
at low knee frequency but should not be used for high knee frequencies as such. 

Our noise bias estimates for EE and BB spectra are equal, but the averaged spectra for the Monte Carlo maps appear visually 
different in this respect. To ensure that this is only due to Monte Carlo noise we ran Madam in Monte Carlo mode, simulating noise 
on the fly and avoiding the costs associated with storing of the time ordered data. After averaging over 117 A^e = 8 noise map 
spectra we found that the differences between EE and BB noise map spectra were at most 10%. 

In the Appendix lA.il we replot some of these figures after subtracting the analytical bias and dividing by the Monte Carlo sample 
deviation to highlight the differences between the analytical and numerical results. 

6.2.3. By power spectrum estimation 

Our final validation procedure for the noise covariance matrices was to use them in Q estimation. Due to resource constraints this 
exercise was conducted at a lower N^q = 8 resolution. All three map-making codes produced maps from the 25 noise realizations. 
Each realization was paired with an independent realization of the CMB sky and the co-added maps were processed using the Bolpol 
code, an implementation of the QML estimator described in Sect 14.11 The 25 power spectrum estimates for each multipole were 
then averaged over and the Bolpol-determined error bars were accordingly divided by V25. 

The example of the results is shown in Fig. Q3J These estimates were obtained for the case with the low knee frequency, 
/knee = 10mHz, using the conventional destriper, Springtide. We note they do not hint unambiguously at any problem with the 
estimated covariance, even if the^ 2 (strongly) and the noise bias (mildly) tests may indicate otherwise. This is likely in part due to 
a lower sensitivity of the power spectrum test on the one hand and on the other due to the fact that the lower resolution has been 
used in this last case. 

Similar statistically good agreements can also be seen in the case of the higher value of /knee, if the covariance is computed 
using either the optimal or generalized destriping technique with the short baselines of 1.25s. If longer baselines are used, i.e., 60s, 
the estimates of the polarized spectra, both E and B, are visibly discrepant with the assumed inputs. Similar disagreement can be 
seen if the off-diagonal elements of the covariance matrices are neglected. In both of these cases, no particular effect on the total 
intensity spectrum can be noticed in the range of investigated angular scales. We illustrate all these statements in Appendix IA.21 
These observations emphasize the importance of precise estimation of the noise covariance in particular for the polarized power 
spectra. 

6.3. Low-resolution maps 

In the previous sections we have discussed our ability to estimate correctly the properties of the noise present in the low-resolution 
maps. We have demonstrated that this is indeed the case for all considered resolution downgrading strategies and both optimal 
and destriper maps. Though in the latter case a baseline length needs to be carefully chosen depending on the time domain noise 
characteristic. 

In this Section we focus on the low-resolution maps themselves. We will look at them from three different perspectives, evalu- 
ating the level of the map-making artifacts left in the maps, the properties of the sky signal and the level of the noise. 

In Fig. [16] we show the differences of the noise-free low resolution maps computed using different downgrading approaches 
discussed earlier and the input map used for the simulations. The reference low resolution version of the input map has been 
obtained via simple binning of the sky signal directly into low-resolution pixels. 

In the case of the direct method we see clearly the extra power spread all over the sky in all three Stokes parameters. Any other 
proposed approach clearly fares much better leading to a substantial decrease of the level of the observed artifacts. This is quantified 
with the help of the pseudo-spectra in Fig.[T7]and in Table[3j where we have collected the root mean square estimates for the signal 
residual maps. 

The CMB part of the low-resolution maps also depends on the downgrading technique. These attempt to suppress the flight 
(subpixel) power and therefore can potentially affect the CMB angular power spectrum even within the band of interest. Fig. [19] 
shows the full-sky pseudo-Q spectra averaged over 117 CMB realizations, downgraded using inverse noise weighting and a number 
of smoothing kernels. Comparison of the spectra shows that it is hard to attain sub-percent bias even at £ = 2N S id Q . If these were 
estimates from an actual power spectrum estimation code with correlated noise and a sky cut, the smoothed covariances would, 
however, be regularized by adding white noise effectively leading to a considerable uncertainty already at that multipole due to 
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ure 15. Averaged power spectrum estimates over 25 noise and CMB realizations. The noise has a 10 mHz knee frequency. 
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Figure 16. Examples of signal striping. We show the difference between a binned and destriped signal-only maps. Rows correspond 
to direct Madam results, 1.25 s and 60 s, and inverse noise weighted 1.25 s baseline case respectively. 
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Figure 17. Comparison of CMB and stripe power spectra reveals that the striping can significantly bias the EE and BB power 
spectrum estimates. Stripe map spectra are computed from maps shown in Fig.QjJ 

white noise and a sky cut. Methods that produce less than 5% bias for Q estimates at t = 2.5A^ S id e are the Gaussian 10° symmetric 
beam and the apodized step function for = (2 or 2.5N S i^ Q , 3N S ^ Q ). 

In Fig. [18] we show the actual sky signal spectra estimated for the low-resolution, noiseless maps. It can be seen that the 
estimated band powers are not drastically affected with respect to the estimates coming from the binned maps. Nevertheless, the 
case of estimates coming from the direct low-resolution maps show some deviation from the binned case. The bias is most prominent 
in the BB low multipole estimates. We note that since the the signal covariance matrix is ill-conditioned it was regularized by adding 
a small white NCM (<x ^ 1 //K) and each signal map received a noise realization consistent with this white NCM. 

None of the proposed downgrading approaches can yield a noise level better than the direct method. This is because noise- 
weighted downgrading or smoothing both introduce departures from the optimal weighting of the noise present in the data. The 
expected level of the noise is therefore an important metric with which to compare the different downgraded maps. Fig. [20l shows 
the analytical noise biases evaluated from an Af s id e =8 Madam noise covariance matrix after smoothing using the 6 different beam 
window functions defined in Sect. 15.41 and the unsmoothed case that is a close match to the noise weighted maps. The narrowest 
Gaussian window function having an FWHM equal to 5°, slightly less than average pixel width, appears pathological due to aliasing 
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ure 18. Band power estimates of CMB with CMB and foreground stripes. 
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Figure 19. Bandwidth limiting the signal using various window functions. In the QML method the quadratic map function is 
multiplied with the inverse Fisher matrix to produce the QML power spectrum estimate (see Sect. I4.lt . The inverse Fisher matrix 
can correct some of the aliasing effects that cause bias in the power spectra. For this figure the pseudo-Q spectra were computed 
from the full sky CMB maps. To simulate the effect of the inverse Fisher in QML we deconvolved our pseudo spectra with a 
mode coupling kernel that we computed from a map of ones. The three panels show the same curves first as raw estimates, then 
after deconvolving the smoothing and pixel windows and finally after subtracting and dividing by input model. The apodized step 
window functions ("cosine") correspond to choices of the thresholds as (20,24), (16,24) and (16,20) respectively. Here, 

solid lines are for Gaussian windows and dashed lines for the apodized step functions. Note that the noise weighting and direct 
low-resolution map-making produce similar aliasing effects. 



effects. The rest of the test cases are more stable but feature a non-negligible amount of aliased power for multipole moments beyond 
£ = 3A^ S ide (the Q are not normalized with the conventional £{l + l)/27r). 
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Figure 20. TT noise bias computed from smoothed covariance matrices. Solid lines correspond to the Gaussian window functions 
and the dashed ones to apodized step functions. Left: Linear vertical scale. Right: Logarithmic vertical scale. 



6.4. Resource requirements 

Table H] lists some CPU time costs for various Af s ide, baseline length and knee frequency combinations. B eing a conside rably ex- 
pensive operation, we have not tested the scaling of the optimal calculation for this particular exercise. See Borrill (1999) for more 
discussion of scaling. The Madam resource cost scales roughly linearly with respect to the number of pixels and the so-called 
baseline correlation length. 

Resource requirements of the three approaches vary. As with map-making, the destriping problem size is related to the chosen 
baseline offset length. The same consideration applies also for noise covariance estimation. Both the Madam generalized destriper 
and the MADping covariance calculations scale by the length of the noise filter. ROMA does the calculation in Fourier space and as 
a results scales as the the logarithm of the noise filter length. The MADping algorithm has no dependency (ignoring communication 
and final file writing) on the number of pixels being used. Although the computation prefactors are not specified, we can see that at 
^V S ide = 32, MADping is already more resource efficient than ROMA (cf. rows two and three). 



28 R. Keskitalo et al.: Residual noise covariance for Planck low-resolution data analysis 

Table 4. CPU time costs for 12 70GHz detector years = 3 • 10 10 samples. Calculations are done on 2.6GHz Quad-core Opterons. 



Case 


side 




PEs 


CPUh 


MADping 


32 


131,073 


14,000 


80,000^ 


MADping 


32 


4,097 


14,000 


7,525 


ROMA 


32 


4,097 


512 


25,000 


ROMA 


4 


16,385 


256 


480 


ROMA 


4 


4,097 


256 


410 


Madam 1.25s, 5(£ 


64 


1,759 


1,024 


901 


Madam 1.25s, 5(T 


32 


1,759 


512 


114 


Madam 1.25s, 10 


32 


8,285 


512 


438 


Madam 10s, 50 


32 


289 


128 


42 


Madam 60s, 50 


32 


84 


64 


34 


Madam 1.25s, 50 


16 


1,759 


256 


34 


Springtide, 10 


32 


1 


1,024 


256 


Springtide, 10 


16 


1 


1,024 


51 



a The filter or baseline correlation length 

b performed on an earlier, dual-core version of the machine 

c 10 and 50 refer to \/f knee frequencies in mHz 



7. Conclusion 

We have presented the formalism and tools to compute the residual noise covariance matrix for three map-making paradigms studied 
for Planck (an optimal method and two destriping methods). The structure of these matrices follows from the scanning strategy but 
is modulated by the underlying noise model that defines the map-making method. The matrices were tested against Monte Carlo 
noise maps that were processed from correlated noise streams into maps using MADmap, Madam and Springtide map-making 
codes. 

The most accurate correspondence between the covariance matrix and the noise maps is, as expected, between the optimal 
map-makers, MADmap and ROMA, and their covariance matrices and the two codes produce nearly identical matrices. Both the 
generalized (Madam) and classical (Springtide) destripers are shown to disregard some medium frequency correlated noise that 
cannot be modelled by the chosen baseline offset length. It is shown that for a low knee frequency, 10 mHz, the Springtide baseline 
length of 1 hour is sufficient to model the correlated noise and compute the residual noise covariance. For a high knee frequency, 
50 mHz, even the Madam 60 s baselines are too long to suffice. However, using a short 1.25 s baselines (just 96 samples) the Madam 
results are extremely close to optimal results even for the high knee frequency. 

As a concluding test we used the matrices in actual power spectrum estimation and verified that all methods model residual 
noise adequately when the noise approximation (baseline length) is short enough to model the correlated noise. 

Resource costs of the methods vary greatly. Although both MADping and ROMA arrive at the same result, the implementations 
differ and the ROMA result scales with the resolution of the map. Both optimal implementations are extremely resource intensive. 
The Madam method can be used to produce good approximations of the optimal covariance matrices at a fraction of their cost. The 
covariance matrices for these tests were evaluated for two low resolutions, N^q = 8 and N^q = 32. It is possible to compute the 
matrices up to N^q = 64 (already 162 gigabytes) or even up to N^q = 128 (2.6 TB) but the computational scaling of the methods 
using the matrices will likely set limits to the usefulness of such resolutions. 

We studied two classes of downgrading strategies, those that make an attempt to limit the signal bandwidth and those that do 
not. The choice of the best downgrading approach depends on both the accuracy of the resulting noise and signal models. First 
measuring our ability to compute an accurate noise covariance matrix and the second describing our ability to control signal effects 
such as striping and aliasing. 

All methods to produce low-resolution maps have their drawbacks. Direct map-making at low resolution produces an unaccept- 
able level of signal striping that is caused by subpixel structure. Downgrading by noise weighting biases the power spectrum through 
aliasing effects. The frequently used Gaussian beam smoothing has a significant drawback of suppressing the signal at otherwise 
useful angular resolutions. We find that an apodized step function is able to retain a great deal of signal power up to £ max = 2Af S id e but 
even then the power spectrum estimates will be biased beyond i - 2.5Af S id e . However, to accurately evaluate the noise covariance 
matrix for a smoothed map, we would need to compute an usmoothed covariance matrix at the high map resolution and then apply 
the same smoothing kernel to both the matrix and the map. Disregarding this requirement leads to disagreement between the map 
and the matrix that can be alleviated by combining two or more of the downgrading methods. 

Of the downgrading methods considered, we consider smoothing, with a suitable choice of the window function and possibly 
a intermediate downgrading step by inverse noise weighting, to produce the best possible low-resolution maps for power spectrum 
analysis. 

We presented in this work a method to compute the residual noise covariance of a smoothed, bandwidth-limited map. The 
method was shown to produce an accurate description of the noise in the smoothed maps when both the map and the matrix agree 
prior to smoothing. Our method of smoothing the covariance matrix makes it possible to consider bandwidth limited low-resolution 
maps and produce sub-percent level unbiased power spectrum estimates up to £ = 2.5Af S id e . 

In this work we have assumed a single frequency channel, uncorrelated noise between detectors, noise that is white at high 
frequencies and full sky coverage. In further work any of these constraints can be lifted. The only application that we used the 
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covariance matrix was in power spectrum estimation. The downgrading methods that suit power spectrum analyses best may not be 
optimal for different low-resolution analysis, e.g., study of large scale topology. A relevant future direction to explore is the use of 
the covariance matrices as inputs in a likelihood code for cosmological parameters. 
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Appendix A: Additional material 

A.1. Noise biases 

For completeness, we present in Fig. IA. li the noise bias computed from the MADmap NCM and the 25 corresponding noise maps 
and in Fig. lA.2l the Madam bias for /k nee = 10 mHz and 60 s baseline. 

To highlight differences between the estimates and simulated noise maps we also show the fractional differences in Figs. IA.3l - 
IA.5I These plots complete the ones presented in Sect. 16.2.21 

A.2. Power spectra 

Here we present another successful test of the covariance matrix used in power spectrum estimation, Fig. IA.6I We also show how 
the power spectrum estimates can be used to pick out inaccurate residual noise co variances in Fig. lA.7[ 
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Figure A.l. Analytical and mean Monte Carlo noise biases from MADmap runs at /knee = 50mHz. Grey band is the l-cr region for 
the average, computed by dividing the sample variance by V25. Like the TE, TB and EB biases are both consistent with zero and 
are not shown here. The analytical bias corresponds to the unsmoothed case presented in Sect.f 
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Figure A.2. Analytical and mean Monte Carlo noise biases from Madam runs at /knee = 10 mHz. 
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Figure A.3. Averaged noise biases after subtracting the analytical estimate and normalizing with the standard deviation. This plot 
contains the same curves as Fig. IA.ll 
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Figure A.4. Averaged noise biases after subtracting the analytical estimate and normalizing with the standard deviation. This plot 
contains the same curves as Fig.[T2l 
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Figure A.5. Averaged noise biases after subtracting the analytical estimate and normalizing with the standard deviation. This plot 
contains the same curves as Fig.[T4l 
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Figure A.6. Averaged power spectrum estimates over 25 noise and CMB realizations. The noise has a 50 mHz knee frequency. 
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Figure A.7. Averaged power spectrum estimates over 25 noise and CMB realizations. The noise has a 50 mHz knee frequency. 



